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Abstract — An efficient, joint transmission delay and channel 
parameter estimation algorithm is proposed for uplinlc asynchronous 
direct-sequence code-division multiple access (DS-CDMA) systems based 
on the space-alternating generalized expectation maximization (SAGE) 
framework. The marginal likelihood of the unknown parameters, 
averaged over the data sequence, as well as the expectation and 
maximization steps of the SAGE algorithm are derived analytically. 
To implement the proposed algorithm, a Markov Chain Monte Carlo 
(MCMC) technique, called Gibbs sampling, is employed to compute the 
a posteriori probabilities of data symbols in a computationally efficient 
way. Computer simulations show that the proposed algorithm has 
excellent estimation performance. This so-called MCMC-SAGE receiver 
is guaranteed to converge in likelihood. 

Index Terms- Asynchronous DS-CDMA, space-alternating generaUzed 
expectation maximization(SAGE), Markov Chain Monte Carlo (MCMC), 
Gibbs sampling. 

I. Introduction 

The performance of direct-sequence code-division multiple-access 
(DS-CDMA) transmission over mobile fading channels depends 
strongly on the reliability of channel parameter and quality of syn- 
chronization for each user: state-of-the-art detection algorithms that 
exploit multiple-access-interference and inter-symbol-interference re- 
quire very powerful estimation algorithms. 

Substantial amount of relevant references appeared in the literature 
on delay estimation. Namely, a new prospective is presented in [1] for 
the maximum likelihood (ML) time-delay estimation. Code timing 
estimation in a near-far environment for DS-CDMA systems was 
introduced in [2]. Joint symbol detection, time-delay and channel 
parameter estimation problems for asynchronous DS-CDMA systems 
have been investigated in several previous works (e.g., [3], [4]). Most 
of these works either work on one signal at a time and treat the other 
signals as interference, or employ a training sequence to obtain a 
coarse estimate of the channel parameters which is consequently used 
to detect data. It is clear that these approaches have disadvantages of 
having higher overhead and additional noise enhancement. 

Some other proposed approaches for joint blind multiuser detection 
and channel estimation for DS-CDMA systems are subspace-based 
and linear prediction-based methods. Subspace-based method usually 
require singular value decomposition or eigenvalue decomposition 
which is computationally costly and does not tolerate mismatched 
channel parameters. Another drawback of this approach is that 
accurate rank determination may be difficult in a noisy environment 
[5], [6]. Moreover, it is not clear how these methods can be extended 
to include the estimation of the transmission delays jointly with the 
channel parameters. 
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The expectation maximization (EM) and space alternating EM 
(SAGE) algorithms are ideally suited to these kind of problems as 
they are guaranteed to converge in likelihood. Earlier work related 
with delay estimation based on the EM algorithm has appeared 
[7], [8]. Efficient iterative receiver structures are presented in [9], 
[10], performing joint multiuser detection and channel estimation 
for synchronous as well as asynchronous coded DS-CDMA systems 
operating over quasi-static flat Rayleigh fading channels, under the 
assumption that the transmissions delays are known. The Bayesian 
EM/SAGE algorithm can be used for joint soft-data estimation and 
channel estimation but the computational complexity of the resulting 
receiver architecture is non-polynomial in the number of users [11]. 
To overcome this draw-back, Hu et al. applied the Variational 
Bayesian EM/SAGE algorithms to joint estimation of the distributions 
for channel coefficients, noise variance, and information symbols 
for synchronous DS-CDMA in [12]. Our work may be considered 
to be a twofold extension of the work by Gallo et al. in [11]: 
First, the proposed receiver performs joint channel coefficient and 
transmission delay estimation within the SAGE framework. Secondly, 
the implication of the Monte-Carlo method in the SAGE framework 
makes it possible to compute soft-data estimates for all users at 
polynomial computational complexity, as well. Here, an efficient 
Markov chain Monte Carlo (MCMC) technique [13] called Gibbs 
sampling is used to compute the a posteriori probabilities (APP) 
of data symbols [14]. The APP's can be computed exactly with the 
MCMC algorithm, which is significantly less complex than a standard 
hidden Markov model approach. The resulting receiver architecture 
works in principal fully blind and is guaranteed to converge. For 
uncoded transmission, a few pilot bits must be inserted, though, to 
resolve the phase ambiguity problem. 

The theoretical framework for the joint transmission delays and 
channel estimation as well as the data detection algorithms can easily 
be extended to coded transmission. 

II. System Description 

We consider an asynchronous single-rate DS-CDMA system with 
K active users using binary phase shift keying (BPSK) modulation 
sharing the same propagation channel. The signal transmitted by each 
user experiences flat Rayleigh fading, which is assumed to be constant 
over the observation frame of L data symbols. Each user employs 
a random signature waveform for transmitting symbols of duration 
Tfc, such that each symbol consists of Nc chips with duration Tc = 
Tt/Nc where Nc is an integer The received signal is the noisy sum 
of all user's contribution, delayed by the propagation delays Tk € 
[0,ri,/2), where the subscript k denotes the label of the fcth user 
After down-converting the received signal to baseband and passing 
it through an integrate-and-dump filter with integration time Ts — 
Tc/Q, Q £ Q^, QNc{L+l) samples over an observation frame of L 
symbols are stacked into a signal column vector r G c<?^<:(^+i)-i 
Note that sampling is chip-synchronous without knowledge of the 



individual transmission delays. It can thierefore be expressed as 

r = S{T)Ad + w. (1) 

In tliis expression the matrix S{t) G £Qr^c(L+i)-ixLK ^.g^i^^^^^ 
the signature sequences of all the users 

S{t) = [ Si(ri), S2{r2), ■■■ ,SK(rK) ] 
where 5fc(rfc) G c'^n,_(l+i)-ixl j^^^ ^j^^ ^^^^ 



Sk{Tk) = 



5fc(rfe,0) 5fe(rfe,l) 



Sk{Tk,L-l) 



and the spreading code vector Sk{Tk,£) G C'^^'^^^"'"^^ ^'*Ms given 
by 



Sk{rkJ) 



Sk{Tk,t) 



The vector Sk{rk,t) contains the spreading code of user k having 
support [lNcT,,{l+l)N,T,] with energy sl{rk,t)Sk{rk,t) = 1- 
Finally, Oa/xi denotes the M x 1-dim. all-zero column vector. 

The block diagonal channel matrix A G C^^'^^^' in (1) is 
given by A = diag{Ai, ■ ■ ■ , Ak}- The channel matrix for user fc, 
Ak G C^^^, is given by Ak = \l® a,k where II is the L-dim. 
identity matrix, and the symbol (d denotes the Kronecker product. The 
fcth user's channel coefficient ak is a circularly symmetric complex 



The maximization (M)-step computes a value of the argument 
in (|2j to obtain the update sj!^^'. The objective function is non- 
decreasing at each iteration. 



B. The Monte-Carlo SAGE algorithm 

We will see that direct computation of the expectation in l|2j 
requires a non-polynomial number of operations in the number 
of users K and thus becomes prohibitive with increasing K. To 
make the computation of the expectation in l|2j feasible though, 
we propose to use the technique of Markov chain Monte Carlo 
(MCMC) to obtain the Monte-Carlo SAGE algorithm. MCMC is a 
statistical technique that allows generation of ergodic pseudo-random 
samples d'*'^' , . . . , d'*''^*' from the current approximation to the 
conditional pdf p(d|r, These samples are used to approximate 
the expectation in l|2j by the sample-mean. The Gibbs sampler and the 
Metropolis-Hastings algorithm are widely used MCMC algorithms. 
Here we describe only the Gibbs sampler [16], [14], as it is the most 
commonly used in applications. Having initialized d'"'"' randomly, 
the Gibbs sampler iterates the following loop at SAGE iteration i: 

• Draw sample dj*'*' from p(di [dj ^\ 



. Draw sample d'''*' fromp(d2|d^''*l , d'''*"'' . . . , df^''"'' , r, 6>W) 
• Draw sample d^'*' from p(dif IdJ*'*' , . . . , d^'^j^, r, fl''') 



Following this approach, we have 

Aft 



user's transmission delay is assumed to be uniformly distributed. Qk{dk,d^^^) = — — |logP fr,d'^'' 

The symbol vector d G f"^^'^ faVp^thp fnrm = rnUW, ... ' ^ ' 



takes the form d = col{di , ■ ■ ■ , (Ik} 



,ak,Tk,a^ 



where the vector dk G C contains the fcth user's symbols, i.e. dk = 
col{dfe(0), ■ ■ ■ ,dk{L - 1)} with dk{l) G denoting the 

symbol transmitted by the fcth user during the ^th signalling interval. 
Finally, the column vector w G ([;'3^c(^+i)-i contains complex, 
circularly symmetric white Gaussian noise having covariance matrix 
A'^ol. We assume that the vectors a = col{ai, 02, ■ ■ ■ t = 

col{ri, T2, ■ ■ ■ , tk}, d and w and their components are independent. 
The receiver does not know the data sequences, the (complex) channel 
coefficients, or the transmission delays. 

III. Monte-Carlo sage Joint Parameter Estimation 

A. The SAGE Algorithm 

In previous applications, the SAGE algorithm [15] has been 
extensively used to iteratively approximate the ML/MAP estimate 
of a parameter vector with respect to the observed data r. To 
obtain a receiver architecture that iterates between soft-data and 
channel estimation, one might choose the parameter vector as 6 — 
{$H(ai), • ■ ■ ,$K(aK), 3(ai), ■ ■ ■ , 3(aA'), n, ■ ■ ■ , tk}. The symbols 
SH(-) and denote the real and imaginary parts of the complex 
argument, respectively. At iteration i, only the parameter vector of 
user fc, 9 k are updated, while the parameter vectors of the other users 
= 0\0k are kept fixed. In the SAGE framework r is referred to 
as the incomplete data. The so-called admissible hidden data Xk ^i^h 
respect to 6 is selected to be Xk ~ {'"'■ '^}- Notice that Xk ^^'^ °nly be 
partially observed. Applying the SAGE algorithm to MAP parameter 
estimation, yields the expectation (E)-step 



Notice that with increasing Nt, the Monte-Carlo SAGE algorithm 
converges to the MAP solution 6 = 6* to random fluctuations 
around 6* [17]. 

C. Receiver design 

This subsection is devoted to the derivation of a receiver architec- 
ture for joint estimation of parameters within the Monte-Carlo SAGE 
framework. Discarding terms independent of a and r, we obtain 

logp(r, d, a, r) = logp(r |d, a, T)+logp(d)+logp(a)+logp(r). 

(3) 

From lO, it follows that 

logp(r|a, T, d) a lR{r'^ SAd} - i/x(6», d) V(f , d), (4) 

where ti{9,d) = Ylk^i^^^o ^k{i,Tk)akdk{t) and (.)''" is the 
conjugate transpose of the argument. 

1 ) The E-step: Substituting Q into ([3j yields after some algebraic 
manipulations for the E-step of the Monte-Carlo SAGE algorithm: 



L-l 



iVc 



Qfc(6'fe,6>W) = Sd{logp(r,d,afc 



1^0 

with the branch definition 



No 



(5) 



(2) 



<f{£,rk) ^ SU£,Tk) (d^:^{£)r ~ll:\e) 



and the interference term 

+Sk' (I - 1, rW) (dfe Wdl7(^ - 1)) " 



Moreover, 



and 



(6) 



xP(dfc(^) = m,dfc,(/) = n I r.T^.aW), for fc' / fc, (7) 

where S = { — 1,+!} is the signal constellation and the lag is within 
range £' € {£-!,£,£ + 1}. 

2) The M-step: The M-step of the SAGE algorithm is realized by 
first maximizing jsjl with respect to the transmission delays Tk, 



Xt, ^ 



(8) 



Then by inserting ([8} into taking derivatives with respect to the 
ak's, setting the results equal to zero, and solving yields 



Ai+i) 



L + No/al 



TV. Monte-Carlo Implementation to the Computation of 
A Posteriori Probabilities 

A. Computation of the soft-data symbols in ([fij 

Let = d\{dk{£)}. For notational simplicity we use d = 

dk{£) throughout this section. Then, the a posteriori probability of 
dk{£) in © can be evaluated as 



P{dk{t)=m\ r,TW,aW) 

= Y^P{dk{£) = m\ d,r,rW,aW) P(d|r,rW,aW) 



1 ^' 

LY,P{d,{£) = m0^''\r,T^^,a^^-^ 



Nt 



(9) 



To compute P{dk{£) = m|d'''*', r, r''' , a''') for this Markov chain 
Rao-Blackwellization technique, we define 



P(rffe(^) = -l|d'''",r,TM,aM)' 

an, the data symbols are i.i.d. a 
)ws from ilOl that 

^l.,t] ^ P(r |dfc(£)^+l,d'''",rW,aW) 
P(r I dk{£) = -l,d''''',TH,aW)' 
from which it can be easily seen that 



(10) 



For uncoded transmission, the data symbols are i.i.d. and equally 
likely. Therefore, it follows from ilOl that 



(11) 



From (TJ, we have p{r\D) ~ exp{ — -^\r — Gd\^), with G 



S{t)A and d — col{di, d2, ■ ■ ■ , dx},- After some algebra Jllb can 
be expressed as 



(12) 



where q = kL + £, and Gq is G with its gth column removed. 
Similarly, dq denotes the vector d with its gth component removed. 

In summary, for each k = 1, 2, • ■ ■ ,K and ^ = 0, 1, • ■ ■ , L — 1, 
to estimate the a posteriori probabilities P(dfc(^) |r, r'*' , a'"') in 
the Gibbs sampler runs over all symbols Nt times to generate a 

collection of vectors which are used in M2\ 

to estimate the desired quantities. 

B. Computation of the soft-value for the product of two data symbols 
in 

=[i,t] \i,t] 
Similarly, a number of random samples d = dk^k'i^') ,t = 

1, 2, • ■ ■ ,Nt,£' £ { — 1, 0, +1} are drawn, using the Gibbs sampling 

technique, from the joint conditional posterior distribution, P(d 

r.T^.aW). Based on the samples d''"'', (dk{£)dk'{£')'^ in Q 
can be evaluated by 

(dki7)d^i£')Y^ ~ (l/Nt) 



Nt 



X E] E/ tnnP ( dk{£) = m,dy{£') = n \ d 

t — 1 m.n^S 



= [i,t] 



We need to evaluate the probability in the expression above. Follow- 
ing the same route taken as in the previous section and after some 
algebra, it can be expressed as 



P{dk{£)^m,dk>{£')^n\d''' 



1 + exp (-<''■•'!) 1 +exp(-A[''«l) 
The quantities C'*'*' and A'"'*' (13} are given by 



(13) 



AMU^K{m(gW)t(r--G!;l 4'")}. 



where p = k'L + £' and q = kL + £. Gp:q is G with its pth and gth 
columns QpjQq removed. Similarly, dp^q denotes the vector d with 
its pth and gth components removed. 

V. Performance ANALYSIS 

A. Modified Cramer-Rao Bounds for the Estimated Parameters 

We now derive the modified Cramer-Rao lower bounds (MCRB) 
on the variances of any unbiased estimates 6 of the parameter vector 
e. It is shown in [18] that for Op G 0, vai(ep - Op) > [7"^(0)]pp, 
where I (6) is the 3K x 3A' Fisher information matrix whose {p, g)th 
component is defined by 

lnp(r, a \ r) 



dOpdOq 



for p, q = 1, 2, • ■ ■ , 3K. 



For the joint likelihood function in (HJl, it is shown in [19] that the 
Fisher information matrix can be computed by 

'dp,\e,d) dn{e,d)' 



(14) 



P,"? = 1,2, • 



,3K. 



Taking the expectations with respect to channel coefficients a and 
data d after taking the partial derivatives in J14t with respect to 6p and 
6q, for different regions of p and q values, under the assumption that 
the data sequences are independent and equally likely and the fact that 
S\Tp,t)S[Tp,t) = 1, for p= 1,2--- /C; ^ = 0, 1,-- - ,L- 1, the 
Fisher information matrix becomes a diagonal matrix whose (p, p)th 
component can be evaluated as 



with the short-cut S' [(] 




S'{1) 



P = K + 1,-- 



,2K 

2K + ,3K. 

(15) 

, . The final result for 



the MCRBs on the estimates of the channel coefficients and the 
transmission delays is obtained by inverting the diagonal matrix 1(0) 
in l ll5b as follows. 



var(afc) > No/L, 
var(?i) > l/(87r^L> 



(16) 
(17) 



= 1, 2, . . . , iC. The symbol 7^ = al/No is the average SNR, B^^ 
is the Gabor bandwidth of the fcth user's spreading code waveform, 
Sfe(i) i.e.. 



B.. = ( / f\S,{f)fdf 



1/2 



and Sk{f) is the Fourier transform of Sk{t),t £ [0,T(,]. Note that 
the Gabor bandwidth Ba^. tends to infinity for rectangular- shaped 
(continuous-time) chip waveforms. 

B. Numerical Examples 

To assess the performance of the proposed (non-linear) Monte- 
Carlo SAGE scheme, an asynchronous uncoded DS-CDMA system 
with K = 5 users, rectangular chip waveforms with processing gain 
Nc = 8, and L = 80 transmitted symbols per block is considered. 
The receiver processes Q = 12 samples per chip. For each data 
block, Gibbs sampling is performed over 50 iterations. A few, say 
Lp — 4 pilot symbols are embedded in each block to overcome 
the phase ambiguity problem. Each user's strongest path from the 
MMSE estimate of a given the K Q {Lp + 1) — 1 samples of r 



The 



and the pilot symbols, yield the initial estimates a'"' and t'"'. 
MMSE estimate of d, given r and weighted by a'"' yields the initial 
symbol estimate d'"'. We refer to this method as MMSE-separate 
estimation (MMSE-SE). For comparison purpose, the SAGE-scheme 
for joint data detection and channel estimation in [10] for known 
transmission delays and hard-decision decoding has also been 
considered subsequently. We refer to these scheme as "SAGE-JDE, 
T known". 

To study the behavior of the proposed MCMC-SAGE scheme, we 
consider communication over AWGN (not known to the receiver). 
The individual powers are given by 

af = -4 dB, cr| = -2 dB, cr| = dB, 
o-| = +2 dB, cr| = +4 dB, 

Fig. [T] shows the mean-square-error (MSE) of the channel estimates 
Si (weakest user) and 03 (normal user) as a function of the 
normalized transmission delays t /Tb which are uniformly distributed 
on the interval between zero and the value on the abscissa. It can be 
seen that the MCMC-SAGE performs close to the MCRB over the 
entire range of t. Not shown in the plot, convergence is achieved after 
around 25 iterations i.e., every user's parameter vector is updated five 
times. 
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Fig. 1. var(afe) of the MCMC-SAGE in near-far scenario. 
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Fig. 2. var(rfc) of the MCMC-SAGE in near-far scenario. 



Fig. [2] depicts the MSE of the delay estimates ri and T3. Notice 
that the MCRB for r tends to zero for time-continuous signature 
waveforms. It can be seen that user 3 does not encounter delay 
estimation errors for small transmission delays i.e., r/Tj, < 0.2. This 
effect can be partially explained by the large number of samples per 
chip i.e., Q = 12. Though for higher transmission delays, var(T3) is 
finite, because of the increasing residual interference in the receiver. 

The bit-error-rate (BER) of the proposed receiver is plotted in 
Fig. [3] versus the effective SNR ^^^k, % = al/No, k = 
1,...,K. The transmission delays are uniformly distributed on 
[0, Tt/2). It can be seen that the MMSE-SDE scheme cannot handle 
delay estimation errors at all due to high correlations between the 
users' signature sequences. The proposed MCMC-SAGE scheme and 
the "SAGE-JDE, r known" perform similar. The weakest user 1 
performs close to the single-user (SU) bound. The normal user 3 
has a multiuser efficiency of roughly 1 dB over the entire range of 
SNR values. 
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Fig. 3. BER-performance in near-far scenario. 
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VI. Conclusions 
A computationally efficient estimation algorithm has been pro- 
posed for estimating the transmission delays and the channel coef- 
ficients jointly in a non-data-aided fashion via the SAGE algorithm. 
The a posteriori probabilities needed to implement the SAGE algo- 
rithm have been computed by means of the Gibbs sampling technique. 
Exact analytical expression have been obtained for the estimates of 
transmission delays and channel coefficients. At each iteration the 
likelihood function is non-decreasing. 
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